This is the second of two Required Practice activities for today. This one focuses on introducing the tmap package. This package is much easier to use than plotly in my view, but it is more specialized. The tmap package is designed for visualizing geospatial data (i.e., mapping) rather than for generating a wide variety of exploratory data analysis graphics (e.g., barcharts, scatterplots).


This intro uses the same shapefile of (redacted) flood insurance claims paid in the state of Virginia compiled by the Federal Emergency Management Administration included in the .zip file. You can read more about the source of these data here if you like. We saw this dataset last week as part of the exam review.


Below, we read in the shapefile from hard disk using the sf::st_read function…


va_nfip_claims_0 <- st_read("./NFIP_Claims_VA.shp")
## Reading layer `NFIP_Claims_VA' from data source 
##   `C:\Users\bw6xs\Documents\PLAN 6122\2023\Lectures\Week 9\Week 9 - Lecture 1 - Required Practice\NFIP_Claims_VA.shp' 
##   using driver `ESRI Shapefile'
## Simple feature collection with 2959 features and 38 fields
## Geometry type: POINT
## Dimension:     XY
## Bounding box:  xmin: -82.8 ymin: 36.5 xmax: -75.3 ymax: 39.4
## Geodetic CRS:  WGS 84
# Filter out records where year of loss is 2003
va_nfip_claims_since_2003 <- va_nfip_claims_0 %>%
  filter(yrOfLss >= 2003)


Now let’s practice making static AND interactive maps with the tmap package 👎. 

# Get a layer with boundaries for context...
bound_0 <- st_read("https://opendata.arcgis.com/datasets/e3c8822a4adc4fc1a542a233893a46d4_0.geojson")
## Reading layer `SDE_USDC_CENSUS_VA_COUNTY' from data source 
##   `https://opendata.arcgis.com/datasets/e3c8822a4adc4fc1a542a233893a46d4_0.geojson' 
##   using driver `GeoJSON'
## Simple feature collection with 134 features and 20 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: -83.67539 ymin: 36.54076 xmax: -75.16644 ymax: 39.46601
## Geodetic CRS:  WGS 84
bound <- st_transform(bound_0, st_crs(va_nfip_claims_since_2003))

# Group and count by county (locality)
by_county <- va_nfip_claims_since_2003 %>%
  drop_na(amnPOBC) %>%
  group_by(contyCd) %>%
  count()


# Join to existing sf object for plotting...
# by_county_joined <- st_join(bound, by_county, left = TRUE)


# Repair faulty feature geometries in the boundary layer...
st_is_valid(bound, reason = TRUE)
##   [1] "Valid Geometry"                               
##   [2] "Valid Geometry"                               
##   [3] "Valid Geometry"                               
##   [4] "Valid Geometry"                               
##   [5] "Valid Geometry"                               
##   [6] "Valid Geometry"                               
##   [7] "Valid Geometry"                               
##   [8] "Valid Geometry"                               
##   [9] "Valid Geometry"                               
##  [10] "Valid Geometry"                               
##  [11] "Valid Geometry"                               
##  [12] "Valid Geometry"                               
##  [13] "Valid Geometry"                               
##  [14] "Valid Geometry"                               
##  [15] "Valid Geometry"                               
##  [16] "Valid Geometry"                               
##  [17] "Valid Geometry"                               
##  [18] "Valid Geometry"                               
##  [19] "Valid Geometry"                               
##  [20] "Valid Geometry"                               
##  [21] "Valid Geometry"                               
##  [22] "Valid Geometry"                               
##  [23] "Valid Geometry"                               
##  [24] "Valid Geometry"                               
##  [25] "Valid Geometry"                               
##  [26] "Valid Geometry"                               
##  [27] "Valid Geometry"                               
##  [28] "Valid Geometry"                               
##  [29] "Valid Geometry"                               
##  [30] "Valid Geometry"                               
##  [31] "Valid Geometry"                               
##  [32] "Valid Geometry"                               
##  [33] "Valid Geometry"                               
##  [34] "Valid Geometry"                               
##  [35] "Edge 2377 has duplicate vertex with edge 2403"
##  [36] "Valid Geometry"                               
##  [37] "Valid Geometry"                               
##  [38] "Valid Geometry"                               
##  [39] "Valid Geometry"                               
##  [40] "Valid Geometry"                               
##  [41] "Valid Geometry"                               
##  [42] "Valid Geometry"                               
##  [43] "Valid Geometry"                               
##  [44] "Valid Geometry"                               
##  [45] "Valid Geometry"                               
##  [46] "Valid Geometry"                               
##  [47] "Valid Geometry"                               
##  [48] "Valid Geometry"                               
##  [49] "Valid Geometry"                               
##  [50] "Valid Geometry"                               
##  [51] "Valid Geometry"                               
##  [52] "Valid Geometry"                               
##  [53] "Valid Geometry"                               
##  [54] "Valid Geometry"                               
##  [55] "Valid Geometry"                               
##  [56] "Valid Geometry"                               
##  [57] "Valid Geometry"                               
##  [58] "Valid Geometry"                               
##  [59] "Valid Geometry"                               
##  [60] "Valid Geometry"                               
##  [61] "Valid Geometry"                               
##  [62] "Valid Geometry"                               
##  [63] "Valid Geometry"                               
##  [64] "Valid Geometry"                               
##  [65] "Valid Geometry"                               
##  [66] "Valid Geometry"                               
##  [67] "Valid Geometry"                               
##  [68] "Valid Geometry"                               
##  [69] "Valid Geometry"                               
##  [70] "Edge 12 has duplicate vertex with edge 835"   
##  [71] "Valid Geometry"                               
##  [72] "Valid Geometry"                               
##  [73] "Valid Geometry"                               
##  [74] "Valid Geometry"                               
##  [75] "Valid Geometry"                               
##  [76] "Valid Geometry"                               
##  [77] "Valid Geometry"                               
##  [78] "Valid Geometry"                               
##  [79] "Valid Geometry"                               
##  [80] "Valid Geometry"                               
##  [81] "Valid Geometry"                               
##  [82] "Valid Geometry"                               
##  [83] "Valid Geometry"                               
##  [84] "Valid Geometry"                               
##  [85] "Valid Geometry"                               
##  [86] "Valid Geometry"                               
##  [87] "Valid Geometry"                               
##  [88] "Valid Geometry"                               
##  [89] "Valid Geometry"                               
##  [90] "Valid Geometry"                               
##  [91] "Valid Geometry"                               
##  [92] "Valid Geometry"                               
##  [93] "Valid Geometry"                               
##  [94] "Valid Geometry"                               
##  [95] "Valid Geometry"                               
##  [96] "Valid Geometry"                               
##  [97] "Valid Geometry"                               
##  [98] "Valid Geometry"                               
##  [99] "Valid Geometry"                               
## [100] "Valid Geometry"                               
## [101] "Valid Geometry"                               
## [102] "Valid Geometry"                               
## [103] "Valid Geometry"                               
## [104] "Valid Geometry"                               
## [105] "Valid Geometry"                               
## [106] "Valid Geometry"                               
## [107] "Valid Geometry"                               
## [108] "Valid Geometry"                               
## [109] "Valid Geometry"                               
## [110] "Valid Geometry"                               
## [111] "Valid Geometry"                               
## [112] "Valid Geometry"                               
## [113] "Valid Geometry"                               
## [114] "Valid Geometry"                               
## [115] "Valid Geometry"                               
## [116] "Valid Geometry"                               
## [117] "Valid Geometry"                               
## [118] "Valid Geometry"                               
## [119] "Valid Geometry"                               
## [120] "Valid Geometry"                               
## [121] "Valid Geometry"                               
## [122] "Valid Geometry"                               
## [123] "Valid Geometry"                               
## [124] "Valid Geometry"                               
## [125] "Valid Geometry"                               
## [126] "Valid Geometry"                               
## [127] "Valid Geometry"                               
## [128] "Valid Geometry"                               
## [129] "Valid Geometry"                               
## [130] "Valid Geometry"                               
## [131] "Valid Geometry"                               
## [132] "Valid Geometry"                               
## [133] "Valid Geometry"                               
## [134] "Valid Geometry"
boundary <- st_make_valid(bound)


# Try again with repaired sf object...
by_county_joined <- st_join(boundary, by_county, left = TRUE)


# Use the qtm function to peruse a basic map WITHOUT having to write much code...
qtm(by_county_joined, fill = "n")

# Set the mode to plot (i.e., static maps) and polish the map a bit...
tmap_mode("plot")

tm_shape(by_county_joined) + 
  tm_polygons("n", palette = "RdYlBu", 
              colorNA = "white",
              textNA = "No claims documented", 
              title = "Number") +
  tm_layout(title = "Flood Insurance Claims Since 2003", frame = FALSE)


This is fine… I guess… 😪 but the strength of tmap is its ability to quickly spin up interactive maps! The first two lines of code in the chunk below display the basemap options that are available for interactive maps generated with the tmap package.


I’ve offered my elderly neighbor $20 to let me try out her stair lift.

I think she’s going to take me up on it. 😆


styles <- schema()$layout$layoutAttributes$mapbox$style$values
styles
##  [1] "basic"             "streets"           "outdoors"         
##  [4] "light"             "dark"              "satellite"        
##  [7] "satellite-streets" "carto-darkmatter"  "carto-positron"   
## [10] "open-street-map"   "stamen-terrain"    "stamen-toner"     
## [13] "stamen-watercolor" "white-bg"
tmap_mode("view")

tm_shape(by_county_joined, name = "NFIP Claims") + 
  tm_polygons(col = "n", palette = "RdYlBu", 
              colorNA = "white",
              textNA = "No claims documented", 
              title = "Number", 
              id = "NAME") +
  tm_layout(title = "Flood Insurance Claims Since 2003", frame = FALSE) +
  tm_basemap(c("Esri.WorldGrayCanvas", "CartoDB.Positron", "Esri.WorldTopoMap"))


Note that by default, the pop-up that appears when we hover over or click on a polygon is the value we assigned to the col argument, but we can format these any way we like… now click a polygon instead of hovering over it 🐭 to see the difference.


tmap_mode("view") 

tm_shape(by_county_joined, name = "NFIP Claims") + 
  tm_polygons(col = "n", palette = "RdYlBu", 
              colorNA = "white",
              textNA = "No claims documented", 
              title = "Number",
              popup.vars = c("Locality name: " = "NAMELSAD", "NFIP claims: " = "n")) +
  tm_layout(title = "Flood Insurance Claims Since 2003", frame = FALSE) + 
  tm_view(set.view = c(-77.4360, 37.5407, 7)) + 
  tm_basemap(c("Esri.WorldGrayCanvas", "CartoDB.Positron", "Esri.WorldTopoMap", "CartoDB.DarkMatter"))
# We can also reverse the palette...
tmap_mode("view")

tm_shape(by_county_joined, name = "NFIP Claims") + 
  tm_polygons(col = "n", palette = "-RdYlBu", 
              colorNA = "white",
              textNA = "No claims documented", 
              title = "Number",
              popup.vars = c("Locality name: " = "NAMELSAD", "NFIP claims: " = "n")) +
  tm_layout(title = "Flood Insurance Claims Since 2003", frame = FALSE) +
  tm_view(set.view = c(-77.4360, 37.5407, 9)) +
  tm_basemap(c("Esri.WorldGrayCanvas", "CartoDB.Positron", "Esri.WorldTopoMap", "CartoDB.DarkMatter"))

If we want to access cool basemaps, it is best to sign up for MapBox token by visiting this page.


Your Turn


Insert a new code chunk that generates at least one static and at least one interactive map of your choice using tmap. This can be with the dataset we have here OR another dataset of your choosing.

When you are done, knit your R Notebook to HTML, then publish the results to RPubs by following the guidelines below:



The above image shows where the Publish button is located in the Viewer pane. Click it and choose Publish HTML.



The next image shows that we are pushing this to RPubs and not to another platform.





Then, we sign in and give the file a name and a description as well as choosing the final bit of the URL.



Finally, if you need to remove something from RPubs, log in and use the Delete button shown here.

Note: you may want to sign up for an account at the outset here.



You have reached the end! 💪